Negative specific heat in a thermodynamic model of multifragmentation 
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We consider a soluble model of multifragmentation which is similar in spirit to many models which 
have been used to fit intermediate energy heavy ion collision data. In this model c v is always positive 
but for finite nuclei c p can be negative for some temperatures and pressures. Furthermore, negative 
values of c p can be obtained in canonical treatment. One does not need to use the microcanonical 
ensemble. Negative values for c p can persist for systems as large as 200 paticles but this depends 
upon parameters used in the model calculation. As expected, negative specific heats are absent in 
the thermodynamic limit. 
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I. INTRODUCTION 



This paper deals with specific heats of an assembly of interacting nucleons. In recent times the subject has 
received a great deal of attention [1-6]. The topic is beset with many controversies. Some of the ideas are : under 
suitable condtions, nuclear systems exhibit negative heat capacities; negative heat capacities are obtainable only in 
the microcanonical ensemble; negative heat capacities also appear in canonical models but disappear once the drop 
■ size crosses the value « 60. 

We investigate the specific heats using a thermodynamic model. The basic assumption of the model is that 
f^i ■ populations of different channels are dictated solely by phase space considerations. This is a common theme in many 
applications, for example, SMM (statistical multifragmentation model) [1] and MMMC (microcanonical metropolis 
Monte Carlo) model [2] although details vary from one model to another. 

A canonical model based on this assumption was shown to be easily soluble requiring only very quick and simple 
computing [7]. The first application used one kind of particle but was later extended to two kinds of particles [8,9]. 
This appears to be accurate enough for many applications [10] and will undoubtedly be used more and more in the 
"^T* ! future. We investigate the question of specific heat in this model primarily using one kind of particle. Two kinds of 
particles were also used but requires longer computing time but we expect no changes from the lessons learnt from 
the model of one kind of particles. We will however show also some results obtained from using two kinds of particles. 

What we will show is that although c v for this model is always positive, c p can sometimes be negative. This is a 
finite particle number effect and negative values disappear in the thermodynamic limit. Furthermore we get negative 
. , values of c p in the canonical model itself. We did not need to go to the microcanonical description. Thermodynamic 
limit is obtained by using the grand canonical ensemble whereas finite systems are described by canonical model 
J-j \ with exact particle number. We find that negative values of c p can persist for fairly large systems although this is 
- - i dependent upon binding energies used etc.. This was not investigated in detail. 

The statements made above appear to hold for the Lattice Gas Model as well. It was demonstrated that c v is 
positve in the Lattice Gas Model even for a very small system. This can be shown almost analytically without having 
to use a Monte-Carlo simulation [4] On the other hand c p is much harder to calculate in the Lattice Gas Model. 
Chomaz et al find that c p can be negative in the Lattics Gas Model [3]. 

For completeness we describe the canonical thermodynamic model in Section II. In the next section we set up the 
grand canonical model to get to the thermodynamic limit. Subsequent sections will show the results. 



II. THE THERMODYNAMIC MODEL 



The thermodynamic model has been described in many places [7,8,11]. For completeness and to enumerate the 
parameters we provide some details. We describe the model for one kind of particles only. The generalisation to two 
kinds could be found in [8,11]. 

If there are A identical particles of only one kind in an enclosure at temperature T, the partition function of the 
system can be written as 
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Qa = jyM- 4 (2.1) 

Here u> is the partition function of one particle. For a spinless particle this is u = ■p-(27rmT) 3 / 2 ; m is the mass of 
the particle; V is the available volume within which each particle moves; A\ corrects for Gibb's paradox. If there are 
many species, the generalisation is 

i 

Here u>i is the partition function of a composite which has i nucleons. For a dimcr i = 2, for a trimer i = 3 etc. 
Equation (2.2) is no longer trivial to calculate. The trouble is with the sum in the right hand side of cq. (2.2). The 
sum is restrictive. We need to consider only those partitions of the number A which satisfy A — irii. The number of 
partitions which satisfies the sum is enormous when A is large. We can call a given allowed partition to be a channel. 
The probablity of the occurrence of a given channel Pin) = P(ni,n 2 , n 3 ....) is 

The average number of composites of i nucleons is easily seen from the above equation to be 

<ni> =^^±± (2.4) 
Qa 



Since i n i = A one readily arrives at a recursion relation [12] 

k=A 

A 



^ k=A 

Qa = -j kuj kQA-k (2.5) 
fc=l 

For one kind of particle, Qa above is easily evaluated on a computer for A as large as 3000 in matter of seconds. It 
is this recursion relation that makes the computation so easy in the model. Of course, once one has the partition 
function all relevant thermodynamic quantities can be computed. 

We now need an expression for uj k which can mimic the nuclear physics situation. We take 

tu k = ^{2vmTf/ 2 k3l 2 q k (2.6) 

where the first part arises from the centre of mass motion of the composite which has k nucleons and q k is the internal 
partition function. For k = 1, q k = 1 and for k > 2 it is taken to be 

q k = exp[(M/ fc - a(T)k 2 ' z + T 2 k/e )/T] (2.7) 

Here, as in [1], Wo=16 MeV is the volume energy term, <r(T) is a temperature dependent surface tension term and 
the last term arises from summing over excited states in the Fermi-gas model. The value of eo is taken to be 16 MeV. 
The explicit expression for a(T) used here is a(T) = (r a [(T 2 - T 2 )/{T 2 + T 2 )] 5 / 4 with <t =18 MeV and T c = 18 
MeV. In the nuclear case one might be tempted to interpret V of eq.(2.6) as simply the freeze-out volume but it is 
clearly less than that; V is the volume available to the particles for the centre of mass motion. Assume that the only 
interaction between clusters is that they can not overlap one another. This assumption restricts the validity of the 
model to low density limit as was stressed in all previous applications of of the model. In the Van der Waals spirit we 
take V — V f reeze — V ex where V ex is taken here to be constant and equal to Vo = A/p . The precise value of V ex is 
inconsequential so long as it is taken to be constant. Calculations employ V; the value V ex enters only if results are 
plotted against p/po = Vq/(V + V ex ) where p is the freeze-out density. 

In the past, calculations with one kind of particle used the parametrisation of eq. (2.7) for all k's however large. 
This means that if the system has A nucleons, the largest possible cluster allowed in the system also has A nucleons. 
While we will show a few cases with this specification we will also consider a variation. We will take the value of q k 
to be given by eq. (2.7) up to a limit k = N and zero afterwards. When A > N, this simply means that the largest 
cluster has N nucleons. 

Using standard definitions: E = T 2 dl g® A and pressure p = T 9L q^ a we arrive at 

E = J2<n k > [E k k m + El nt ) (2.8) 
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where Ef n = §T and E l k nt = k(-W + T 2 /e ) + cr(T)fc 2 / 3 - T[da(T)/dT}k 2 / 3 . The last term in E l k nt was neglected 
in [7]. It is included here but makes little difference. The expression for pressure is 

P=^J2<n i > (2.9) 

Multiplicity m is given by m = < n i >■ 

III. INFINITE MATTER LIMIT: GRAND CANONICAL MODEL 

If there were only monomers in the grand canonical ensemble one would solve 

p = expfa/Tfa (3.1) 

where (Dj = u>i/V of section II. Given p and T one then finds the chemical potential p. The number of particles is then 
given by A — pV where V and A are very large (thermodynamic limit). The fluctuations in the number of particles 
implied by the use of grand canonical ensemble are then negligible compared to the average number A. 

If we have a model where the only allowed species are monomers and dimers and the total particle number is very 
large one would solve: 

p = exp( / u/T)o;i + 2 x cxp(2yu/T)a)2 (3.2) 

where phase-space consideration has implied that chemical equilibration exists, that is, the chemical potential of the 
dimer is twice that of the monomer, i.e., p 2 = 2/z. 

For a system which is very large but, for which, the heaviest cluster has TV nucleons and no more, one needs to 
solve 

N 

P = ^2kexp(kfj,/T)Q k (3.3) 
i 

In this case, one might argue that one is considering a model in which the composites obey eq.(2.6) up to k = N and 
Wfe's for k > N arc all zeroes. Of course it is possible that both A and N are very large. Use of the grand canonical 
ensemble always implies that A is very large but N may be large or small. 

Pressure in the grand canonical model is calculated from p = (T/V)lnZ granc [ which in this model reduces to 
p = {T/V)^2 < n i > where < > /V = exp(ifj,/T)Qi. Notice that formally this eq. is same as eq.(2.9) but, of 
course, < rij > in cq. (2.9) is calculated according to the canonical formula, eq.(2.4). 

IV. SPECIFIC HEATS IN THE MODEL 

In [7] where the canonical thermodynamic model was first studied for phase transitions, it was pointed out that for 
a given density p, the specific heat per particle c„ = Cy/A tends to oo at a particular temperature when the particle 
number A tends to oo. Since, 

r —(<>E\ — T ( §H\ — T ( 9 2 f \ 

— \ OTlv ~ \dT>V ~ \WF ) ' y ' 

a singularity in Cy signifies a break in the first derivative of F, the free energy and a first order phase transition. 
The specific heat c v in the model has been studied in more than one application and always found to be positive. 
We now turn our attention to c p studied in this canonical model. It is instructive to look at p — p curves at different 
temperatures (equation of state) to gain an understanding (Fig.l). For 200 particles (A=200 and 7V=number of 
nucleons of the largest allowed cluster=200) this is drawn at three temperatures: T\ < T 2 < T 3 . Here T 2 is only 
slightly higher than T\. We notice that on T\, T 2 isothermals there are regions of mechanical instability where dp/dp 
is negative. It is in this region that one encounters negative values for c p . Instead of p let us use the variable V cx 1/p. 
Thus regions of mechanical instability are characterised by dp/dV >0. Let us try to understand how this can happen. 
In simpler cases as in a gas of noninteracting monomers, the multiplicity m (which determines the pressure, i.e., 
p = T x y) is simply A and (dp/dV)T is always negative, (we actually use m — 1 rather than m for calculating p but 
this is immaterial for our discussion). In the thermodynamic model, because of composites, m << A at moderate 
temperatures. At fixed temperature T > 0, m will always increase with V. Negative compressibility is marked by 
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(w) T > V- -^ et us consider points a, b (region of positive c p ) and points c, d (region of negative c p ). Points b, c are 
on I\ = T and points a, d are on T 2 = T + <5T with > 0. Using 
p = Tf = (T + (5T)f±|f we arrive at 

5m <5V ST 

In the region (c, d), (5V is negative, ST ia always positive thus 5m is negative. If m goes down then so does the kinetic 
energy and also the potential energy (creating more m creates more surface and hence more energy). Thus in this 
region with increasing temperature but constant pressure, both the kinetic and the potential energy of the system 
go down. In the "normal" region Sm is positive and both the kinetic and the potential energies increase with T at 
constant pressure. This is illustrared in Table 1. Finally we show, in fig. 2, the caloric curve for a given pressure 
where in part of the curve temperature does go down with excitation energy. The fall is very gentle whereas the rise 
with energy when it happens is faster. 

The occurrence of a negative C p in spite of a positive Cy is allowed in the following well-known relation [17]: 

C p ~ C v = VT^ 

where a is the volume coefficient expansion and k is the isothermal compressibility given by: 

1 fdV 

a= v{df /r 

For negative k, C p is less than Cy and can become negative. 
Using the equality, 

(dv\ -_(dv\ (dp\ 
\dT) p - \dp) T \dT) p 

we can also write, 

C p -Cv = VT v ± (§£) p . 

This shows that, C p can drop below Cy if isobaric volume coefficient of expansion becomes negative which is the case 
in some regions of Fig.l. 



V. EXTRAPOLATION TO THERMODYNAMIC LIMIT 



For A very large but N=200 we use the grand canonical ensemble. For a given p and T, we solve for p where N in 
eq. (3.3) is set at 200. In the p — p diagram there are no regions of mechanical instability (see fig. 3) . For comparison, 
the p — p diagram for N — 200 but ^4=200 obtained by the canonical calculation is also shown in the same figure. We 
see that in the low density side (the gas phase) the two diagrams coincide. The rise of pressure with density is quite 
rapid and linear. After the two diagrams separate, the rise of pressure with density in the grand canonical model 
slows down considerably but there is no region of mechanical instability although the canonical calculation with 200 
particles has a region of instability. In the grand canonical result which represents the thermodynamic extrapolation, 
we have not reached the classic liquid-gas coexistence limit where there would be no rise of pressure at all (like in 
Maxwell construction). We think the reason is this. The largest cluster has size 200 which is not a big enough 
number. Condensation into the largest and larger clusters still does not behave like a liquid. We now increase the 
largest cluster size to 2000. Now the coexistence region is very clear and there is unmistakable signature of first order 
phase transition. This grand canonical result is very close to the case where N = oo ( see Fig. 3 in [13]). In the same 
figure we also show results of a canonical calculation with A=2000 and ^=2000. The region of mechanical instability 
has gone down considerably but it has not disappeared showing that we have not reached the thermodynamic limit 
yet. 
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VI. MORE CALCULATIONS IN THE CANONICAL MODEL 



The mechanical instability which led to negative values of c p is not only a finite number effect but it is also 
dependent on details of parameters; see also [6]. As in fig. 3 we draw a p — p diagram for 200 particles in fig. 4 but 
now the largest cluster has N = 100, that is, ivk is given by eq. (2.6) up to fc=100 and is zero for k > 100. The 
mechanical instability region has completely disappeared. In fact, the negative compressility in the p — p diagram 
in fig. 4 disappears even with the following minimal change. We use qk of cq. (2.7) upto k — 100 and for k > 100 
use qk = exp[0.97(Wofc — er(T)fc 2 / 3 + T 2 fc/e )/T]. We were surprised that with such small changes zones of negative 
compressibility disappeared but such consequences were anticipated in other models before [6] . 



VII. CHEMICAL POTENTIALS 



In this section we deal with two kinds of particles and discuss the behaviours of chemical potentials as the proton 
fraction of large systems changes. This is remotely connected with specific heats but the behaviour of chemical 
potentials with the proton fraction has attracted some attention in recent times and we felt that it is of general 
interest to show what the behaviour is in the themodynamic model. It was shown that in mean-field theories of 
nuclear matter there is chemical instability as a function of y = Z j (N + Z) in limited regions of y, that is, (dp,p/dy) p ^T 
becomes negative in some region of \ip — y plane (correspondingly (cfytjv /dy) Pt T becomes positive). This is analgous to 
mechanical instability as a function of density [14-16]. We have seen that in the thermodynamic model there arc no 
regions of mechanical instability for large systems (the grand canonical results). We will see that there is no chemical 
instability either in the model in the large particle number limit. Now we need to consider the thermodynamic model 
for two kinds of particles. For details we refer to [8]. A composite has two indices: i=proton number, j=neutron 
number with a = i + j. Analogous to eq.(2.4) we have 

Qz-i,N-j 1 v 

< >= - (7-1) 

HZ,N 

where the nuclear properties are contained in u>ij: 

w id = ^(2nmT)^ 2 a^ 2 qiJ (7.2) 
We take the internal partition function of the composite to be 

qij = cxp[(Wa - era 2 / 3 - s ^—J^- + aT 2 /e )/T] (7.3) 

As is usual in all infinite matter case calculations, the coulomb interaction is switched off. We take W — 15.8 MeV, 
(7 = 18 MeV,s=23.5 MeV and eo=16.0 MeV. For a > 5 we use this formula. For lower masses we simulate the no 
coulomb case by setting the binding energy of 3 He=binding energy of 3 H and binding energy of 4 Li=binding energy 
of 4 H. 

For a given a what are the limits on i(or j = a — i)7 This is a non-trivial question. In the results we will show, we 
have taken the limits by calculating the drip lines of protons and neutrons as given by the binding energy formula. 
Limiting oneself to within the drip lines is a well-defined prescription, but is likely to be an underestimation since 
resonances show up in particle-particle correlation experiments. On the other hand, for a given a, taking limits of i 
from to a is definitely an overestimation. 

In fig. 5 we have drawn isothermals (at T=6.0 MeV) for two component nuclear systems for different y's in the 
grand canonical ensemble. We restrict y between 0.3 and 0.5 and p/ po between and 0.5, the ranges for which the 
thermodynamic model is expected to be reliable. In the calculation the largest cluster is taken to be a = 500. The 
same figure also shows the behaviours of yUp and p,N at constant pressure. The derivative is seen to be always 
positive (simultaneously is negative). In the next figure, for completeness, we have continued the model beyond 
p/po=0.5 and gone upto the highest possible limit of p/po = 1 in the model to see the behaviour of pp and pn- No 
chemical instability is seen. 
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VIII. SUMMARY AND DISCUSSION 



We have shown that with usual concepts one can obtain a negative value of Cp in part of the T — E plane within 
the framework of a thermodynamic model. Although we have shown this, for the sake of simplicity, using one kind 
of particle only, we have checked that the phenomenon remains when a more complicated version with two kinds of 
particles and realistic binding energies for the composites are used. The Cy is positive and its origin is the cost in 
surface energy to break large clusters into smaller clusters and nucleons. A negative C p is seen in our exactly soluble 
canonical ensemble model for small systems. This negative value arises in regions of mechanical instability where the 
isothermal compressiblity is negative or equivalently, the isobaric volume expansion coefficient is negative. A negative 
isobaric volume expansion leads to a decrease in multiplicity, or total number of clusters, with temperature and a 
corresponding decrease in energy For larger systems these regions disappear and in the grand canonical limit, C p is 
always positive. 

Since several papers have demonstrated the existence of negative specific heats it is pertinent to mention the 
relevance of our work to these earlier works. Our model is not, in any simple way, connected to negative specific heat 
found in ten dimensional Potts Model [2]. The specific heat considered in that work is Cy and the negative specific 
heat appears only in microcanonical treatment. The negative specific heat seen here is at least partially similar to 
that seen in [3,4]. There is no negative specific heat in Cy [4] but it makes its presence felt when Cp is considered 
[3]. Although our model is quite different from the one considered in [6] the results are similar. 

Unfortunately we can not recommend any experiments to verify the conclusions of this paper. Nuclear disasembly in 
heavy ion collisions can not be fine tuned. There is no reason to think that it takes place exactly at constant volume or 
exactly at constant pressure. Calculations at constant volume gives quite reasonable predictions for observables that 
have been measured [10,18] but this does not rule out the possibility that variations happen. If disassembly always 
took place at constant pressure then the following idealised experiment would be useful. One measures excitation 
energy per particle and also the temperature. One would then find there are cases where the average excitation energy 
per particle goes down even though the temperature rises. 
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TABLE I. 


Variation of energies per particle (MeV) with temperat 


ure {MeV) 


in the negative and positive 


compressibility 


zones, for p 


= 0.017 MeV /m~ 3 












T 


/9 / On 










6.0 


0.146 


0.978 


-5.235 


-4.257 




6.1 


0.212 


0.638 


-6.970 


-6.332 


6.2 


0.392 


0.294 


-8.708 


-8.414 




6.0 


0.104 


1.422 


-3.271 


-1.849 


dp * u 


6.1 


0.090 


1.653 


-2.513 


-0.859 


6.2 


0.082 


1.824 


-2.027 


-0.202 
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FIG. 1. EOS in the canonical model for a system of A=200. The largest cluster also has iV=200. 
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FIG. 2. Caloric curve at a constant pressure (p — 0.017 MeV fm 3 ) in the canonical model with ^4=200 and iV=200. The 
solid and dashed portions of the curve give -ve and +ve c p respectively. 
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FIG. 3. EOS at T = 6 MeV in the two models. For the left panel the largest cluster has iV=200 and for the right panel 
iV=2000. For the canonical calculation, the left and right panel has A=200 and 2000 respectively, but for the grandcaonocal 
calculations, ^4 = 00. (Sec text). 
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FIG. 4. EOS in the canonical model for a system of 200 particles, but the number of nucleons of the largest cluster is 
restricted to 100. 
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FIG. 5. The top panel are isothermals at T — QMeV for different y's (proton fraction). The lower panels show the behavior 
of hn and fip as a function of y in the density range of the top panel. Calculations are done in a grandcanonical model with 
the largest cluster having 500 nucleons. 



12 








I I I I I I I I I I I I I I I I I I I I I 

0.3 0.35 0.4 0.45 0.5 

proton fraction (y) 



FIG. 6. The behaviors of pn and pp extrapolated to higher densities. For this figure the p — p/po diagram (upper panel 
Fig. 5) was extended upto p/po=l. 
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